home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_12_01 / letters / rpft.c < prev    next >
C/C++ Source or Header  |  1993-06-07  |  22KB  |  674 lines
  1. /******************************************************
  2.  * RPFT.C -- Curve fitting with optional extrapolation.
  3.  * Fits either polynomial or rational function. Based
  4.  * on algorithms defined in Press, et al., "Numerical
  5.  * Recipes", 2nd ed., Cambridge University Press, 1992
  6.  * Developed using the Borland C compiler.
  7.  *
  8.  * To compile:    bcc    rpft.c
  9.  *
  10.  *   Lowell Smith
  11.  *   3368 Crestwood Drive
  12.  *   Salt Lake City, UT  84109-3202
  13.  *
  14.  *  June 8, 1993
  15.  *    1)  Corrected error in COMMD for reading command
  16.  *        line -DIG value
  17.  *    2)  Revised RATINT to eliminate occasional error
  18.  *        interpolating exact polynomial data with a Y
  19.  *        value (not the first or last point) of zero.
  20.  *****************************************************/
  21. #include <stdio.h>
  22. #include <math.h>
  23. #include <string.h>
  24. #include <dir.h>
  25. #include <ctype.h>
  26. #include <stddef.h>
  27. #include <stdlib.h>
  28.  
  29. #define NPFAC 8
  30. #define PI 3.141592653589793
  31. #define PIO2 (PI/2.0)
  32. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  33. void ratlsq(double xs[],double fs[],int npt,int mm,
  34.       int kk, double cof[], double *dev);
  35. double *dvector(long nl, long nh, const char *ner);
  36. double **dmatrix(long nrl,long nrh,long ncl,long nch,
  37.          const char *ner);
  38. void free_dvector(double *v, long ncl);
  39. void free_dmatrix(double **m, long nrl, long ncl);
  40.  
  41. /*****************************************************/
  42. void main(int argc, char **argv)
  43. { double a,b,*xs,*fs,cof[42],dev=0.,sumd,x,yy;
  44.   int i,j,nn,nd,ncof,dig,npt=0,xln,yln;
  45.   char nar[90];
  46.   char filenam[80];
  47.   void commd(int argc, char**argv, double *a,double *b,
  48.         int *nn, int *nd, int *dig, int *xln,
  49.         int *yln, char *filenam);
  50.   void inpt(int nn,int npt,double a, double b, int xln,
  51.         int yln, double *xs,double *fs,char filenam[]);
  52.   void pcshft(double a,double b,double d[],int n);
  53.   void chebpc(double c[],double d[],int n);
  54.  
  55.   commd(argc,argv,&a,&b,&nn,&nd,&dig,&xln,&yln,
  56.         filenam);
  57.   if (xln) { a=log(a);  b=log(b); }
  58.   if (nd)            /* Fit rational function */
  59.   { ncof=nn+nd+1;   npt=NPFAC*ncof;
  60.     xs=dvector(1L,(long)npt,"xs in main");
  61.     fs=dvector(1L,(long)npt,"fs in main");
  62.     inpt(0, npt, a, b, xln, yln, xs, fs, filenam);
  63.     ratlsq(xs,fs,npt,nn,nd,&cof[0],&dev);
  64.     fprintf(stderr,"Est. max error = %.7G\n",dev);
  65.     free_dvector(xs,1L);
  66.     free_dvector(fs,1L);
  67.   }
  68.   else               /* Fit polynomial */
  69.   { ncof=nn+5; xs=dvector(0L,(long)ncof,"xs in main");
  70.     fs=dvector(0L,(long)ncof,"fs in main");
  71.     inpt(ncof, npt, a, b, xln, yln, xs, fs, filenam);
  72.     chebpc(fs,cof,nn+1);
  73.     pcshft(a,b,cof,nn+1);
  74.     free_dvector(xs,0L); free_dvector(fs,0L);
  75.   }
  76.   for (i=0;i<=nn+nd;++i)
  77.   { sprintf(nar,"%.*G ",dig,cof[i]);
  78.     sscanf(nar,"%lf",&cof[i]);
  79.   }
  80.   if (nd) printf("\n      Rational function coefficien"
  81.      "ts\n   Numerator               Denominator\n\n");
  82.   else printf("\n  Polynomial\n  Coefficients\n\n");
  83.   for (i=0;i<max(nd,nn)+1;++i)
  84.   { if (i<=nn) printf("%-24.16G",cof[i]);
  85.     else printf("%*s",24," ");
  86.     if (i<nd) printf("%-24.16G\n",cof[i+nn+1]);
  87.     else printf("\n");
  88.   }
  89.   for(;;)            /* Calculate trial values */
  90.   { i=-1;
  91.     fprintf(stderr,
  92.           "Enter E to quit or a trial X value: ");
  93.     i=scanf("%lf",&x);
  94.     if (i<1) exit(1); if (xln) x=log(x);
  95.     yy=cof[nn]; for (j=nn-1;j>=0;j--) yy=yy*x+cof[j];
  96.     if (nd)
  97.     { for (sumd=0.,j=nn+nd;j>=nn+1;j--)
  98.                  sumd=(sumd+cof[j])*x;
  99.       yy=yy/(1.0+sumd);
  100.     }
  101.     if (yln) yy=exp(yy); if (xln) x=exp(x);
  102.     fprintf(stderr,"For X=%.8G  Y=%.8G\n",x,yy);
  103.   }
  104. }
  105. /******************************************************
  106.  * COMMD - Parses the command line
  107.  *****************************************************/
  108. void commd(int argc, char**argv, double *a,double *b,
  109.         int *nn, int *nd, int *dig, int *xln,
  110.         int *yln, char *filenam)
  111. { struct ffblk ffblk;
  112.   int i,j,k,l;
  113.   void help(char *msg);
  114.  
  115.   *a=1.11e300; *b=-1.11e300; *nn=-32768; *nd=-32768;
  116.   *dig=7, *xln=0, *yln=0;
  117.   if (argc<2) help(""); filenam[0]=0;
  118.   for (i=1,k=0;i<argc;k=0,++i)
  119.   { for (j=0; j<(l=strlen(argv[i])); ++j)
  120.     { argv[i][j]=toupper(argv[i][j]);
  121.       if (argv[i][j]=='=') ++k;
  122.     }
  123.     if (k)
  124.     { if (k>1 || (!isdigit(argv[i][l-1])
  125.                       && argv[i][l-1] !='.'))
  126.       { fprintf(stderr,"Invalid command line parameter"
  127.                  " %s\n",argv[i]);
  128.         help("");
  129.       }
  130.       if (!strncmp(argv[i],"LL=",3))
  131.       { k=sscanf(&argv[i][3],"%lf",a);
  132.         if (k<1) help("Invalid value for command line"
  133.                    " parameter LL\n");  continue;
  134.       }
  135.       if (!strncmp(argv[i],"UL=",3))
  136.       { k=sscanf(&argv[i][3],"%lf",b);
  137.         if (k<1) help("Invalid value for command line"
  138.                    " parameter UL\n");  continue;
  139.       }
  140.       if (!strncmp(argv[i],"ND=",3))
  141.       { k=sscanf(&argv[i][3],"%d",nd);
  142.         if (k<1) help("Invalid value for command line"
  143.                    " parameter ND\n");  continue;
  144.       }
  145.       if (!strncmp(argv[i],"NN=",3))
  146.       { k=sscanf(&argv[i][3],"%d",nn);
  147.         if (k<1) help("Invalid value for command line"
  148.                    " parameter NN\n");  continue;
  149.       }
  150.       if (!strncmp(&argv[i][0],"-DIG=",5))
  151.       { k=sscanf(&argv[i][5],"%d",dig);
  152.         if (k<1) help("Invalid value for optional"
  153.                    " command line parameter DIG\n");
  154.         continue;
  155.       }
  156.       fprintf(stderr,"Invalid command line parameter "
  157.                    "%s\n",argv[i]);
  158.       help("");
  159.     }
  160.     else
  161.     { if (!strncmp(&argv[i][0],"-XLN",4))
  162.       { *xln=1; continue; }
  163.       if (!strncmp(&argv[i][0],"-YLN",4))
  164.       { *yln=1; continue; }
  165.       if (!filenam[0] && argv[i][0] != '-')
  166.       { strcpy(filenam,argv[i]);
  167.         if (findfirst(filenam,&ffblk,0))
  168.         { fprintf(stderr,"Input file %s doesn't "
  169.             "exist\n",filenam); help("");
  170.         }
  171.       }
  172.       else
  173.       { fprintf(stderr,"Invalid command line parameter"
  174.                    " %s\n",argv[i]);   help("");
  175.       }
  176.     }
  177.   }
  178.   k=0;
  179.   if (*a>1.1e300)
  180.   { fprintf(stderr,"Parameter LL undefined\n"); ++k; }
  181.   if (*b<-1.1e300)
  182.   { fprintf(stderr,"Parameter UL undefined\n"); ++k; }
  183.   if (*nn==-32768)
  184.   { fprintf(stderr,"Parameter NN undefined\n"); ++k;
  185.     *nn=5;
  186.   }
  187.   if (*nd==-32768)
  188.   { fprintf(stderr,"Parameter ND undefined\n"); ++k;
  189.     *nd=5;
  190.   }
  191.   if (filenam[0]==0)
  192.     { fprintf(stderr,"Input file is not defined\n");
  193.       ++k;
  194.     }
  195.   if (*nn<1 || (*nd==0 && *nn>20) || (*nd && *nn>10))
  196.   { fprintf(stderr,"Invalid value %d for "
  197.             "command line parameter NN\n",*nn); ++k;
  198.   }
  199.   if (*nd<0 || *nd>10)
  200.   { fprintf(stderr,"Invalid value %d for "
  201.            "command line parameter ND\n",*nd); ++k;
  202.   }
  203.   if (*a>*b)
  204.   { fprintf(stderr,"Lower limit > upper limit\n");
  205.     ++k;
  206.   }
  207.   if (*dig<1 || *dig>20)
  208.   { fprintf(stderr,"Invalid value %d for "
  209.           "command line option DIG\n",*dig); ++k;
  210.   }
  211.   if (*xln && *a<=0.)
  212.   { fprintf(stderr,"Lower limit is <=0  with "
  213.         "logarithmic transform of X values\n"); ++k;
  214.   }
  215.   if (k) help("");
  216. }
  217. /******************************************************
  218.  * HELP - Prints the help screen
  219.  *****************************************************/
  220. void help(char *msg)
  221. { char *hlp[] =
  222.   { "USAGE: rpft LL=a UL=b NN=n ND=m [-DIG=n -XLN -YLN"
  223.     "] file","","   LL=a   a is the lower limit of the"
  224.     "region for fit","   UL=b   b is the upper limit o"
  225.     "f the region for fit","","   ND=m | m>0: rational"
  226.     " function, denominator degree m,","   NN=n |     "
  227.     " numerator degree n.  1<=m<=10  1<=n<=10","      "
  228.     "  | m=0: Polynomial degree n.  1<=n<=20",""," -DI"
  229.     "G=n   (optional) n = number of significant digits"
  230.     ,"              for coefficients on output. Defaul"
  231.     "t=7.","   -XLN   (optional) Transform X values to"
  232.     " ln(X)","   -YLN   (optional) Transform Y values "
  233.     "to ln(Y)","","   file   File with input X Y data "
  234.     "points, 1 per line","","All command line paramete"
  235.     "rs except DIG, XLN and Y